For all of the above-mentioned variational formulations, in order to solve them with the finite element method, it is necessary to generate a system of linear equations. In particular, it is necessary to generate integrals that form a system of linear equations. Remembering that the rows and columns in the global system of equations correspond to the underlying functions - in our example they correspond to the B-spline functions - we generate the following system of equations
\( \begin{bmatrix} & \color{red}{B^x_{1,p} B^y_{1,p}} & \cdots &\color{red}{ B^x_{N_x,p}B^y_{1,p}} \\ \color{blue}{ B^x_{1,p}B^y_{1,p}}& a(\color{blue}{ B^x_{1,p}B^y_{1,p}}, \color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{1,p}B^y_{1,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \color{blue}{ B^x_{2,p}B^y_{1,p}} & a(\color{blue}{ B^x_{2,p}B^y_{1,p}},\color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{2,p}B^y_{1,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \vdots & \vdots & \vdots & \vdots \\ \color{blue}{ B^x_{N_x,p}B^y_{N_y,p} } & a(\color{blue}{ B^x_{N_x,p}B^y_{N_y,p}}, \color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{N_x,p}B^y_{N_y,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y } \\ \end{bmatrix} \)
\( = \begin{bmatrix} l(B^x_1(x)*B^y_1(y)) \\ l(B^x_1(x)*B^y_2(y)) \\ \vdots \\ l(B^x_{N_x}(x)*B^y_{N_y}(y)) \\ \end{bmatrix} \)
B-spline functions that do not have a common support (which are defined over distant parts of the grid) generate zero integrals, and zero values in the matrix. Therefore, when generating a matrix, we should not iterate over the rows and columns of the matrix, but we should iterate over the elements of the computational grid, and over the basis functions that are defined on a given element. The resulting contributions of the integrals should be aggregated into the global matrix. Then our generation of integrals will avoid all non-zero terms of the matrix.
The algorithm contains six nested loops:
1 for \( nex=1,N_x \) (loop through elements along the axis \( x \))
2 for \( ney=1,N_y \) (loop through elements along the axis \( y \))
3 for \( ibx1=1,p+1 \) (loop through \( p+1 \) B-splines of the element along the axis \( x \))
4 for \( iby1=1,p+1 \) ( loop through \( p+1 \) B-spline of the element along the axis \( y \))
5 \( i = f(nex,ibx1) \) ( calculate the B-spline index along the axis \( x \))
6 \( j = f(ney,iby1) \) ( calculate the B-spline index along the axis \( y \))
7 \( irow = g(nex,ibx1,ney,iby1) \) (global index of 2D B-spline)
8 \( l(irow)+= l(B^x_{i,p}B^y_{j,p}) \)
9 for \( ibx2=1,p+1 \) ( loop through \( p+1 \) B-splines of the element along the axis \( x \))
10 for \( iby2=1,p+1 \) (loop through \( p+1 \) B-splines of the element along the axis \( y \))
11 \( k = f(nex,ibx2) \) (calculate the B-spline index along the axis \( x \))
12 \( l = f(ney,iby2) \) (calculate the B-spline index along the axis \( y \))
13 \( irow = g(nex,ibx2,ney,iby2) \) (global index of 2D B-spline)
14 \( M(irow,icol)+= a(B^x_{i,p}B^y_{j,p},B^x_{k,p}B^y_{l,p}) \)
The function \( f(nex,ibx) \) computes the index of the one-dimensional B-spline basis function for element local function \( ibx \) defined along the axis
\( x \) for segment \( nex \). For uniform order B-splines \( p \) this index is equal to \( nex+ibx-1 \).
Functions \( g(nex,ibx,ney,iby) \) calculate the global index of a two-dimensional B-spline. In the case of uniform B-splines this formula is \( irow = i*j \) where \( i=f(nex,ibx) \) and \( j=f(ney,iby) \).
In turn, the computation of integrals is generally performed by numerical integration using points and weights of quadrature, for example Gauss quadrature. The one-dimensional integral of a polynomial on an interval can be calculated exactly as the sum of the values of this polynomial at selected points of the quadrature spread over the interval, multiplied by appropriate weights. The derivation and theory of numerical quadratures is beyond the scope of this manual.
Gauss quadratures are named after a German mathematician who invented this method of integration in 1826 [1].
The choice of quadrature depends on the type of basis function, and it is possible to reduce the number of quadrature points by using special quadratures dedicated to particular basis functions. In general, the basis functions on a single element are polynomials and the use of Gauss quadratures is justified.
The following example illustrates the formation of points and Gaussian quadrature weights on an interval
\( [0,1] \) which guarantees the exact calculation of the integral function of two-dimensional B-splines of a degree not greater than 10.
m_points[1]= (1.0-0.973906528517171720077964)*0.5;
m_points[2]= (1.0-0.8650633666889845107320967)*0.5;
m_points[3]= (1.0-0.6794095682990244062343274)*0.5;
m_points[4]= (1.0-0.4333953941292471907992659)*0.5;
m_points[5]= (1.0-0.1488743389816312108848260)*0.5;
m_points[6]= (1.0+0.1488743389816312108848260)*0.5;
m_points[7]= (1.0+0.4333953941292471907992659)*0.5;
m_points[8]= (1.0+0.6794095682990244062343274)*0.5;
m_points[9]= (1.0+0.8650633666889845107320967)*0.5;
m_points[10]= (1.0+0.9739065285171717200779640)*0.5;
m_weights[1]=0.0666713443086881375935688*0.5;
m_weights[2]=0.1494513491505805931457763*0.5;
m_weights[3]=0.2190863625159820439955349*0.5;
m_weights[4]=0.2692667193099963550912269*0.5;
m_weights[5]=0.2955242247147528701738930*0.5;
m_weights[6]=0.2955242247147528701738930*0.5;
m_weights[7]=0.2692667193099963550912269*0.5;
m_weights[8]=0.2190863625159820439955349*0.5;
m_weights[9]=0.1494513491505805931457763*0.5;
m_weights[10]=0.0666713443086881375935688*0.5;
The quadratures are originally tabulated for the interval \( [-1,1] \) so here we scale them to a range \( [0,1] \) adding 1 to the coordinates of the points and dividing the coordinates of the points by 2 (that is, multiplying by 0.5).
Then, the numerical integration performed on line 14 (and subsequent ones) in the above code requires two additional loops after Gaussian's quadrature points and weights, and it looks like this:
15 for \( m=1,p/2+1 \) (loop through Gaussian quadrature points along the axis \( x \) for B-splines of order \( p \))
16 for \( n=1,p/2+1 \) (loop through Gaussian quadrature points along the axis \( y \) for B-splines of order \( p \))
17 \( M(irow,icol)+=\\ a(B^x_{i,p}(m\_points[m])B^y_{j,p}(m\_points[n]), B^x_{k,p}(m\_points[m])B^y_{l,p}(m\_points[n])) \\ *m\_weights[m]*m\_weights[n]*area \)
In particular, for our example problem of heat transport in an L-shaped area we have
\( M(u,v) = \left(\frac{\partial B^x_{i,p} }{\partial x }(m\_points[m])\frac{\partial B^x_{k,p} }{\partial x }(m\_points[m])\right)*m\_weights[m] \\ + \left(\frac{\partial B^y_{j,p} }{\partial y }(m\_points[n]) \frac{\partial B^y_{l,p} }{\partial y }(m\_points[n]) \right)*m\_weights[n] \)
In the above formula, we calculate the derivative values of the one-dimensional B-spline functions at the points of Gaussian quadrature \( (m\_points[m],m\_points[n]) \) multiplied by the weights of Gaussian quadrature \( *m\_weights[m]*m\_weights[n] \).